Federico Cerisola (cerisola@df.uba.ar)
Departamento de Física, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
Como repasamos en clase, para un oscilador unidimensional amortiguado (con una fuerza de amortiguamiento proporcional a la velocidad) tenemos la siguiente ecuación de movimiento
$$ \ddot{\psi} = -\omega_0^2\psi -\Gamma\dot{\psi} $$con $\omega_0^2 = \frac{k}{m}$ y $\Gamma$ la constante de amortiguamiento por unidad de masa.
Como vimos, la forma funcional de las soluciones cambia según la realción entre $\omega_0$ y $\Gamma$, obteniendo tres casos distintos: subarmotiguado, sobreamortiguado y amortiguamiento crítico.
En este caso, vimos que la solución más general es de la forma
$$ \psi(t) = A\,e^{-\frac{\Gamma}{2}t}\cos(\omega_1 t) + B\,e^{-\frac{\Gamma}{2}t}\sin(\omega_1 t) $$con $\omega_1 = \sqrt{\omega_0^2 - \left(\frac{\Gamma}{2}\right)^2}$. Imponiendo condiciones iniciales vimos además que las constantes $A$ y $B$ son tales que
$$ \psi(t) = e^{-\frac{\Gamma}{2}t}\left( \psi(0)\cos(\omega_1 t) + \left(\dot{\psi}(0) + \frac{\Gamma}{2}\psi(0)\right)\frac{\sin(\omega_1 t)}{\omega_1}\right) $$A continuación se muestra código para graficar las soluciones para disintas condiciones iniciales.
# definimos una funcion que calcula la solucion subamortiguada
# la funcion toma como input el tiempo, omega_0, Gamma
# y las condiciones iniciales en posicion (psi0) y velocidad (vel0)
def psi_subamort(t, omega0, Gamma, psi0, vel0):
# por simplicadad definimos Gamma/2 como una propria variable
Gamma2 = Gamma/2
# definimos el omega1 dado
omega1 = np.sqrt(omega0**2 - Gamma2**2)
psi = np.exp(-(Gamma2)*t)*(psi0*np.cos(omega1*t) + (vel0 + Gamma2*psi0)*np.sin(omega1*t)/omega1)
return psi
# definimos una funcion auxiliar para graficar
# esto nos sirve para despues poder usar la
# interfaz interactiva
def plot_psi_subamort(omega0, Gamma, psi0, vel0):
# vamos a graficar hasta t_max (5 ciclos)
t_max = 5 * 2*np.pi/omega0
# definimos vector con tiempos discretizado en 2000 puntos
t = np.linspace(0, t_max, 2000)
# calculamos solucion
psi = psi_subamort(t, omega0, Gamma, psi0, vel0)
# graficamos usando matplotlib
plt.figure()
plt.plot(t, psi)
plt.xlabel('$t$')
plt.ylabel('$\psi$')
plt.grid()
plt.show()
# creamos widget iteractivo
widgets.interact(plot_psi_subamort,
omega0=widgets.FloatSlider(min=0, max=10, step=1, value=1),
Gamma=widgets.FloatSlider(min=0, max=10, step=1, value=0.5),
psi0=widgets.FloatSlider(min=-1, max=1, step=0.1, value=0.5),
vel0=widgets.FloatSlider(min=-1, max=1, step=0.1, value=0.1))
interactive(children=(FloatSlider(value=1.0, description='omega0', max=10.0, step=1.0), FloatSlider(value=0.5,…
<function __main__.plot_psi_subamort(omega0, Gamma, psi0, vel0)>
En este caso, vimos que la solución más general es de la forma
$$ \psi(t) = A\,e^{-\left(\frac{\Gamma}{2} - \omega_2\right)t} + B\,e^{-\left(\frac{\Gamma}{2} + \omega_2\right)t} $$con $\omega_2 = \sqrt{\left(\frac{\Gamma}{2}\right)^2 - \omega_0^2}$.
Es más, utilizando el truco propuesto en la guía y usando la definición de seno y coseno hiperbólicos, la solución más general, con las condciones iniciales ya impuestas, se puede reescribir como (ej. 2)
$$ \psi(t) = e^{-\frac{\Gamma}{2}t}\left( \psi(0)\cosh(\omega_2 t) + \left(\dot{\psi}(0) + \frac{\Gamma}{2}\psi(0)\right)\frac{\sinh(\omega_2 t)}{\omega_2}\right) $$A continuación se muestra código para graficar las soluciones para disintas condiciones iniciales.
# definimos una funcion que calcula la solucion sobreamortiguada
# la funcion toma como input el tiempo, omega_0, Gamma
# y las condiciones iniciales en posicion (psi0) y velocidad (vel0)
def psi_sobreamort(t, omega0, Gamma, psi0, vel0):
# por simplicadad definimos Gamma/2 como una propria variable
Gamma2 = Gamma/2
# definimos el omega2 dado
omega2 = np.sqrt(Gamma2**2 - omega0**2)
psi = np.exp(-(Gamma2)*t)*(psi0*np.cosh(omega2*t) + (vel0 + Gamma2*psi0)*np.sinh(omega2*t)/omega2)
return psi
# definimos una funcion auxiliar para graficar
# esto nos sirve para despues poder usar la
# interfaz interactiva
def plot_psi_sobreamort(omega0, Gamma, psi0, vel0):
# vamos a graficar hasta t_max (50 veces tiempo 2/Gamma)
t_max = 100 / Gamma
# definimos vector con tiempos discretizado en 2000 puntos
t = np.linspace(0, t_max, 2000)
# calculamos solucion
psi = psi_sobreamort(t, omega0, Gamma, psi0, vel0)
# graficamos usando matplotlib
plt.figure()
plt.plot(t, psi)
plt.xlabel('$t$')
plt.ylabel('$\psi$')
plt.grid()
plt.show()
# creamos widget iteractivo
widgets.interact(plot_psi_sobreamort,
omega0=widgets.FloatSlider(min=0, max=10, step=1, value=0.5),
Gamma=widgets.FloatSlider(min=0, max=10, step=1, value=2),
psi0=widgets.FloatSlider(min=-1, max=1, step=0.1, value=0.5),
vel0=widgets.FloatSlider(min=-1, max=1, step=0.1, value=0.1))
interactive(children=(FloatSlider(value=0.5, description='omega0', max=10.0, step=1.0), FloatSlider(value=2.0,…
<function __main__.plot_psi_sobreamort(omega0, Gamma, psi0, vel0)>
En este caso, vimos que la solución más general, con las condiciones iniciales ya impuestas, es de la forma
$$ \psi(t) = e^{-\frac{\Gamma}{2}}\left[\psi(0) + \left(\dot{\psi}(0) + \frac{\Gamma}{2}\psi(0)\right)t\right] $$A continuación se muestra código para graficar las soluciones para disintas condiciones iniciales.
# definimos una funcion que calcula la solucion de amoritguamiento critico
# la funcion toma como input el tiempo, Gamma
# y las condiciones iniciales en posicion (psi0) y velocidad (vel0)
def psi_amortcrit(t, Gamma, psi0, vel0):
# por simplicadad definimos Gamma/2 como una propria variable
Gamma2 = Gamma/2
psi = np.exp(-(Gamma2)*t)*(psi0 + (vel0 + Gamma2*psi0)*t)
return psi
# definimos una funcion auxiliar para graficar
# esto nos sirve para despues poder usar la
# interfaz interactiva
def plot_psi_amortcrit(Gamma, psi0, vel0):
# vamos a graficar hasta t_max (50 veces tiempo 2/Gamma)
t_max = 100 / Gamma
# definimos vector con tiempos discretizado en 2000 puntos
t = np.linspace(0, t_max, 2000)
# calculamos solucion
psi = psi_amortcrit(t, Gamma, psi0, vel0)
# graficamos usando matplotlib
plt.figure()
plt.plot(t, psi)
plt.xlabel('$t$')
plt.ylabel('$\psi$')
plt.grid()
plt.show()
# creamos widget iteractivo
widgets.interact(plot_psi_amortcrit,
Gamma=widgets.FloatSlider(min=0, max=10, step=1, value=2),
psi0=widgets.FloatSlider(min=-1, max=1, step=0.1, value=0.5),
vel0=widgets.FloatSlider(min=-1, max=1, step=0.1, value=0.1))
interactive(children=(FloatSlider(value=2.0, description='Gamma', max=10.0, step=1.0), FloatSlider(value=0.5, …
<function __main__.plot_psi_amortcrit(Gamma, psi0, vel0)>
Si ahora consideramos además un forazado externo $F(t)$ tenemos la ecuación de movimiento
$$ \ddot{\psi} = -\omega_0^2\psi -\Gamma\dot{\psi} + \frac{F(t)}{m} $$Como en la clase, nos vamos a concentrar en el caso de un forzado armónico $F(t) = F_0\cos(\Omega t)$. Por lo tanto tenemos
$$ \ddot{\psi} = -\omega_0^2\psi -\Gamma\dot{\psi} + \frac{F_0}{m}\cos(\Omega t) $$La solución homogénea es alguna de las soluciones de la sección amortiguada anterior, que cumplen todas que $\psi_{hom}(t) \rightarrow 0$, cuando $t \rightarrow \infty$. Por lo tanto para tiempos largos tenemos que la solucón se parece siempre a la solución particular (régimen estacionario), $\psi(t) \approx \psi_p(t)$.
Para la solución particular mostramos que se puede escribir de la forma
$$ \psi_p(t) = A_{el}\cos(\Omega t) + A_{abs}\sin(\Omega t) $$A continuación se muestra código para graficar la solución particular en función del tiempo.
# función para clacular la amplitud elastica del oscilador forazado
def Ael(m, omega0, Gamma, F0, Omega):
# definimos variables axuiliares
dw = omega0**2 - Omega**2
og = Omega * Gamma
# numerador
num = F0 * dw
# denominador
den = m * (dw**2 + og**2)
return num / den
# función para clacular la amplitud abosrbente del oscilador forazado
def Aabs(m, omega0, Gamma, F0, Omega):
# definimos variables axuiliares
dw = omega0**2 - Omega**2
og = Omega * Gamma
# numerador
num = F0 * og
# denominador
den = m * (dw**2 + og**2)
return num / den
# definimos una funcion auxiliar para graficar
# esto nos sirve para despues poder usar la
# interfaz interactiva
def plot_psi_particular(m, omega0, Gamma, F0, Omega):
# vamos a graficar hasta t_max (5 ciclos de Omega)
t_max = 5 * 2*np.pi/Omega
# definimos vector con tiempos discretizado en 2000 puntos
t = np.linspace(0, t_max, 2000)
# calculamos solucion
psi = Ael(m, omega0, Gamma, F0, Omega)*np.cos(Omega*t) + Aabs(m, omega0, Gamma, F0, Omega)*np.sin(Omega*t)
# graficamos usando matplotlib
plt.figure()
plt.plot(t, psi)
plt.xlabel('$t$')
plt.ylabel('$\psi$')
plt.grid()
plt.show()
# creamos widget iteractivo
widgets.interact(plot_psi_particular,
m=widgets.FloatSlider(min=0.1, max=10, step=0.1, value=1),
omega0=widgets.FloatSlider(min=0, max=10, step=1, value=1),
Gamma=widgets.FloatSlider(min=0, max=10, step=1, value=0.5),
F0=widgets.FloatSlider(min=-1, max=1, step=0.1, value=0.5),
Omega=widgets.FloatSlider(min=0, max=100, step=0.5, value=0.5))
interactive(children=(FloatSlider(value=1.0, description='m', max=10.0, min=0.1), FloatSlider(value=1.0, descr…
<function __main__.plot_psi_particular(m, omega0, Gamma, F0, Omega)>
A continuación se muestra código para graficar los $A_{el}$, $A_{abs}$ y la potencia absorbida en función de $\Omega$.
# función para clacular la amplitud potencia absorbida del oscilador forazado
def Pabs(m, omega0, Gamma, F0, Omega):
Abs = Aabs(m, omega0, Gamma, F0, Omega)
P = F0 * Omega * Abs / 2
return P
# definimos una funcion auxiliar para graficar
# esto nos sirve para despues poder usar la
# interfaz interactiva
def plot_psi_resonancia(m, omega0, Gamma, F0, OmegaMin, OmegaMax):
# definimos vector con Omega discretizado en 2000 puntos
Omega = np.linspace(OmegaMin, OmegaMax, 2000)
# calculamos solucion
El = Ael(m, omega0, Gamma, F0, Omega)
Abs = Aabs(m, omega0, Gamma, F0, Omega)
P = Pabs(m, omega0, Gamma, F0, Omega)
# graficamos usando matplotlib
plt.figure()
plt.plot(Omega, El)
plt.xlabel('$\Omega$')
plt.ylabel('$A_{el}$')
plt.grid()
plt.show()
plt.figure()
plt.plot(Omega, Abs)
plt.xlabel('$\Omega$')
plt.ylabel('$A_{abs}$')
plt.grid()
plt.show()
plt.figure()
plt.plot(Omega, P)
plt.xlabel('$\Omega$')
plt.ylabel('$P_{abs}$')
plt.grid()
plt.show()
# creamos widget iteractivo
widgets.interact(plot_psi_resonancia,
m=widgets.FloatSlider(min=0.1, max=10, step=0.1, value=1),
omega0=widgets.FloatSlider(min=0, max=100, step=1, value=10),
Gamma=widgets.FloatSlider(min=0, max=10, step=0.1, value=0.5),
F0=widgets.FloatSlider(min=-1, max=1, step=0.1, value=0.5),
OmegaMin=widgets.FloatSlider(min=0, max=100, step=0.01, value=0.01),
OmegaMax=widgets.FloatSlider(min=0, max=100, step=0.01, value=20))
interactive(children=(FloatSlider(value=1.0, description='m', max=10.0, min=0.1), FloatSlider(value=10.0, desc…
<function __main__.plot_psi_resonancia(m, omega0, Gamma, F0, OmegaMin, OmegaMax)>